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We demonstrate that an inverse Monte Carlo approach allows to reconstruct effective interaction 
potentials from real-space images. The method is exemplified on monomolecular ethanol-water 
films imaged with scanning force microscopy (SFM), which provides the spatial distribution of the 
molecules. Direct Monte Carlo simulations with the reconstructed potential allow for obtaining 
characteristics of the system which are unavailable in the experiment, such as the heat capacity 
of the monomolecularly thin film, and for a prediction of the critical temperature of the demixing 
transition. 


PACS numbers: 68.15.+e, 05.10.-a, 05.50.+q, 02.50.Tt 

Real-space imaging is widespread in physical, chemi¬ 
cal and biological sciences and is often the only or the 
best way to obtain information about a system. Exam¬ 
ples include, but are not limited to scanning probe or 
optical microscopies. A particular challenge is to predict 
the behavior of an experimental system outside of the 
experimentally accessible set of parameters such as time 
or temperature. This can be achieved given the knowl¬ 
edge of the interactions in the system, which can indeed 
be restored from real-space images, as we discuss in the 
present work. 

This work is motivated by recent investigations of 
a monomolecular film of nanophase separated water- 
ethanol mixtures confined in a slit-pore between a mica 
substrate and a graphene coating as observed by scan¬ 
ning force microscopy (SFM) [TJ [2]. The system ex¬ 
emplifies a situation where the distribution of species 
is known experimentally, yet the interactions between 
them are difficult to estimate. An effective interaction 
emerges from the complex interplay of intermolecular 
forces, graphene’s elastic energy, and (screened) dipole- 
dipole interactions due to charge transfer [2]. In Ref. [2] 
it has only been possible to provide a qualitative dis¬ 
cussion, which has shown that the experimental situa¬ 
tion resembles the behavior of a system with short-range 
attractive and long-range repulsive interactions between 
molecules of the same type, at a temperature above the 
critical point. However, such a discussion could give nei¬ 
ther the strength nor the range of corresponding inter¬ 
actions. Now we use this system as an example for the 
reconstruction of interaction potentials from real-space 
images, which gives the key for a prediction of properties 
which were not immediately observed, including the crit¬ 
ical temperature for the phase separation. The ability to 
estimate the critical temperature is of crucial importance 
for designing new experiments on this interesting system. 


*Electronic address: sokolov@physik.hu-berlin.de 




FIG. 1: (Color online) Binarized SFM topography image: 
White (black) pixels are ethanol-rich (water-rich) areas, re¬ 
spectively. Single-layer graphene, 512 x 512 pixel, 500 nm 
wide. See [2] for more details. The inset (red frame) shows a 
simulated configuration with an effective potential from IMC 
on a 128 x 128 lattice (periodic boundaries). 


The SFM image of a system is a rasterized and pixelled 
image representing the height profile of the graphene film. 
The height profile observed shows a clear bimodal struc¬ 
ture, well-approximated by two Gaussian peaks of similar 
widths [2]; the finite width can be attributed to thermal 
and instrumental noise. The difference in height corre¬ 
sponding to the peaks of the distribution is around 1 A 
and is of the order of the difference in size of water and 
ethanol molecules. The corresponding domains are inter¬ 
preted as water-rich and alcohol-rich clusters. The SFM 
images are therefore filtered and reduced to 1 bit color 
depth (see [2] for details). One of such images is shown 
in Fig. Here the system is in thermodynamic equilib¬ 
rium. 
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The present work is devoted to the question: Which 
quantitative information can be drawn from images as 
the one in Fig.[^ Our approach is based on the determi¬ 
nation of an effective pair interaction potential directly 
from the real-space images by use of the inverse Monte 
Carlo (IMC) method. 

The pixelled SFM images trivially map to square lat¬ 
tices, and hence a lattice binary mixture model suggests 
itself for a coarse-grained description of the situation (the 
coarse graining in experiment is a trivial consequence of 
the relatively low lateral resolution of the SFM, being of 
the order of a few nanometers, so that single molecules 
are not resolved). The model is equivalent to a lattice 
spin system [3] . Such spin models (Ising models with ad¬ 
ditional long-ranged potentials) have been discussed in 
the literature and simulated, see, e. g., [4]. We assign 
the values = +1 = — 1) to the white (black) pix¬ 

els of the image. The translation from the language of 
binary mixture into the one of lattice spin models is dis¬ 
cussed in Appendix IMC is a well-established method 
to calculate effective pair interaction potentials from pair 
correlation functions (CFs), see Ref. [5] for a review, and 
is typically applied in a continuous setting. 

Moreover, we pose the IMC method into a larger scope 
of maximum likelihood estimates. This can allow in fu¬ 
ture for its application to more complex (say, anisotropic) 
situations, to situations under large experimental errors, 
and for the comparison of different models using infor¬ 
mation criteria. For the discussion of the maximal like¬ 
lihood approaches in statistical physics see, e. g.. Ref. [6] 
and references therein. This work is however devoted to 
a different popular inverse problem of statistical physics, 
the inverse Ising model. 

As discussed in [2], the effective interaction between 
the molecules in a slit-pore between the mica surface 
and the graphene film contains essentially three contri¬ 
butions: The direct van-der-Walls interaction between 
the molecules, the effective electrostatic interaction due 
to the charge transfer, and the additional interaction due 
to the elastic deformation of the graphene sheet pushed 
towards the sublayer by the Casimir force and conform¬ 
ing to the molecular relief. The third contributions might 
contain non-additive (multiparticle) components, which 
are however short-ranged (the elastic interaction poten¬ 
tial behaves like 1/r^ 13 i) and decay considerably at 
the distances of the order of the lateral resolution of the 
image. Therefore, in our coarse-grained picture, one can, 
at least as the first approximation, consider only the pair¬ 
wise interactions. The Hamiltonian of a spin system for 
a given spin configuration {cr} is 

^({^}) = E V - ^ext M({a}), (1) 

ijeG 

where W is the effective pair potential, i^ext is an exter¬ 
nal magnetic field and M({cr}) = total 

magnetization of the lattice. Here Vij is the distance 
vector between pixels i, j; G is the set of all pixels/lattice 
sites. 


We moreover introduce the spin-spin correlation func¬ 
tion at some distance vector r on the lattice 

Cr(M) = ^ (2) 

i,3eG\rij=r 

where N = \G\ is the number of pixels/lattice sites. Now 
we can rewrite the Hamiltonian as 


'Hiia}) = NYWr a({<T}) - ifext M(M), (3) 

r 

with Wr = W{r). Eq. § shows that the energy of 
the system depends only on the correlation function C 
and the magnetization M. The case i^ext = 0 corre¬ 
sponds to symmetric black-white coverage, which is (ap¬ 
proximately) exhibited by many - but not all - images. 
The possibility to include i^ext 7^ 0 is discussed in Ap¬ 
pendix 

The probability of a given spin configuration {a} for a 
system with known interaction potential W is given by 
the Gibbs distribution 


p{{a}\W) 


exp [-I3nw{{cr})] 


Zw 


( 4 ) 


where 'Hw{{o'})) is the energy of {a} for given IF, p = 
Ijk^T^ and Zw is the corresponding partition function 
defined as the sum X]{cr} [~{{(^})] over all spin 
configurations. Substituting Eq. (§ into Eq. @ we see 
that the distribution p({cr}| IF) is of the exponential class, 
and that (for i^ext = 0) Cr is a sufficient statistics, i. e., 
gives all information necessary to infer HV, see, e. g., [9]. 

We now consider IF as a vector of parameters (with 
entries Wr)^ and introduce the likelihood £(IF;{cr}) = 
p{{(t}\W). The corresponding log-likelihood L(W) = 
ln£(IF; {a}) is then 


L{W) = -pl-Lw{{(^})-^^Zw = /3[Fw-Hw({c^})], (5) 


where now Fw is the free energy of the spin system for 
given IF at temperature T which is independent of the 
configuration of spins. If we assume that the spin config¬ 
uration observed in the experiment is representative^ the 
energy is well approximated by the internal energy 
Ew of the system, in which case L{W) = ^{Fw — Ew) = 
—Sw/^B, where Sw is the system’s entropy: The log- 
likelihood function for the interaction potential is pro¬ 
portional to the negative value of the entropy of the sys¬ 
tem with a corresponding spin configuration under the 
potential assumed. Its maximization then assumes that 
the most likely value of the interaction potential is such, 
that the experimental configuration of spins provides the 
maximal information about it. We note that the equiva¬ 
lence of the maximum likelihood and minimum entropy 
in statistics is known m- Here, however, the statistical 
entropy has the meaning of the thermodynamical entropy 
pertinent to a representative configuration. 

The representation given by Eq. allows for perform¬ 
ing an IMC procedure using the algorithm proposed in 
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FIG. 2: (Color online) Radial histogram of the “spin-spin” 
CF for the SFM image in Fig. Inset: A part of the CF 
used in the calculation and the CF restored by IMC. 



m, adapted to a lattice system. We consider the CF 
obtained from the experimental image (7, see Fig. as 
a vector of data (with entries (7^), where r has been 
binned radially. The position of the (global) minimum 
at r 9px corresponds to a characteristic size of the 
clusters. We use /cp = 1 in the following. Let us now 
fix some specific distance r and consider the necessary 
condition for the test potential W to maximize L(W) (or 
to minimize Sw) provided the data: 


dL{W) _ dUw 1 dZw 
dWr ^ dWr ^ dWr 


with 


1 dZ\Y 
^ dWr 


W} 


( 7 ) 


Since dl-LwIdWr = N Crdcr}), the sum in the last ex¬ 
pression can be rewritten as the average over the canoni¬ 
cal distribution, i. e., as the canonical correlation function 
Cr{W) = (C'^({cr}))|w^ = T,{a}^r{{cr})p{{cr}\W) for the 
given interaction potential. The brackets (...) and (...)|w 
denote the average over the canonical distribution, the 
latter specifying the dependence on W. We thus see that 
the necessary condition for the extremum of the likeli¬ 
hood is the system of equations 


Cr{W) =Cr for all r. (8) 


If the canonical CF is obtained by the Monte Carlo (MC) 
approach, this inverse problem is solved by IMC. We note 
that fitting the CF by IMC methods is a known approach 
m- However, its validity is typically founded in Hen¬ 
derson’s theorem [12] stating that, provided only pair¬ 
wise interaction, the CF defines the interaction potential 
uniquely, up to the integration constant. 


The necessary condition for the maximum likelihood 
can be rewritten as \/wSw = 0. The algorithmic im¬ 
plementation of the IMC starts from taking some ini¬ 
tial value of the potential vector (close enough to 

the fix point VF* delivering the minimum entropy). For 
W = we want to determine 5W such that 

VwSw = ^swSw = 0 expand this expression to 
the first order in SW. The p-th component of S/swSw 
is given by 


1 

PN dSWp 


= Cp 




+13 N Cov(Cp(M),Cq(M))|,y(o) 

Q 


( 9 ) 


with Cov(^, B)|,y = {AB)\^- {A)\yy The data 

vector Cp is known, the vector Cp(VT*-°^) and the ele- 
ments of the covariation matrix are obtained in a direct 
MC simulation for a given potential Note that the 

emergence of the covariance matrix of Cp stresses again 
the nature of the IMC as a statistical inference method: 
This matrix is proportional to the Fisher information ma¬ 
trix of the inference problem [6]. 

The above system of linear equations ® is solved 
numerically, yielding the (presumably small) correction 
terms SW to the potential. The first iteration is then 
W^^^ = + SW. These updates are repeated un¬ 
til converges to VF*. The method is basically 

a multi-dimensional Newton-Raphson method to solve 
\/wSw = 0, where coefficients are obtained in a sta¬ 
tistical simulation. In our MC simulations we use the 
Glauber dynamics (local spin-flip updates) that violates 
the conservation of particle numbers, in contrast to the 
(slower) Kawasaki dynamics (local particle-exchange up¬ 
dates), which however makes no difference when only 
equilibrium properties are discussed. 

If is already close enough to the solution VF* 

the algorithm converges, but if it is too far, instabili¬ 
ties occur. To enforce convergence for a wider range of 
initial potentials. Ref. [TT] proposed scaling the correc¬ 
tion terms, i. e., using ^ KSW^'^^ with 

A G [0,1] at iteration n + 1. Our experience for the SFM 
images shows that A = 1 almost always leads to instabil¬ 
ities, while suitable values lie in the range A G [0.1,0.5]. 
In order to keep the computational complexity accept¬ 
ably low we considered W to have a rotational symmetry 
(i. e. the same value for all lattice sites at a given sep¬ 
aration; the lattice remains, of course, anisotropic) and 
introduced a cutoff radius Rent that limits the range of 
the pair potential, i. e., Wr = 0 for all r > Rent- 

The iterations are stopped when the difference between 
experimental and simulated correlation becomes small. 
To quantify this, we compute the weighted mean square 
of (7 —C(kF), where the weights are the inverse statistical 
uncertainties in C at the respective distances, and com¬ 
pare it to unity. All details of implementation are given 
in Appendix [B] 

For SFM images with a white-pixel concentration of 
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FIG. 3: (Color online) Reconstructed potential from the cor¬ 
relation in Fig. Parameters: 64 x 64 lattice, Rcut = 22, 
= —0.001. Hand-tuned initial potential, convergence 
after 14 inverse MC iterations with progressively increasing 
measurement times, 7.9 x 10^ global MC updates in total. 


50% (or close to it), F7ext can be fixed to 0 (or some 
suitable small value). For images where this is not the 
case (finite magnetization), we extended the algorithm 
to further reconstruct the correct value of ii^ext- The ex¬ 
tension is in complete analogy to the previous derivation, 
see Appendix [B] 

We tested the procedure on a two-dimensional ferro¬ 
magnetic Ising model with additional long-ranged repul¬ 
sive 1/r^ interaction, as discussed before in the context 
of ultrathin magnetic films [4] which exhibit a clustered 
state at intermediate temperatures. The IMC procedure 
is able to reconstruct the underlying pair-potential from 
CFs obtained from a direct simulation, provided suffi¬ 
cient statistical precision of the input data, see also Ap¬ 
pendix 

The result of the application of the procedure to the 
experimental CF, Fig. is shown in Fig. Inside the 
cutoff radius Rent this potential reproduces the experi¬ 
mental correlation very accurately, cf. the inset in Fig.|^ 
It is attractive at short distances and repulsive at larger 
ones, in agreement with our qualitative expectations. We 
would like to point out that only Cr for r < Rent (range 
of inset in Fig. were used as input data. 

To restore experimentally relevant properties, like the 
critical temperature, we used the reconstructed potential 
in Fig. for direct MC simulations on 64 x 64, 128 x 128 
and 192 x 192 lattices, giving very similar results. We 
examined the specific heat capacity 




- {n? 


( 10 ) 



FIG. 4: (Color online) MC simulations with reconstructed po¬ 
tential (Fig. [^: cn in dependence on T with hxed i4ext = 0. 
Measured with decreasing T in four independent series (dif¬ 
ferent markers) on a 64 x 64 lattice, 10^ MC updates per mea¬ 
surement (of which 10% for thermalization; T G [0.81,0.87] 
was resolved with 10^ MC updates per measurement). 


as well as (specific) magnetization, and magnetic suscep¬ 
tibility. The behavior of cn{T) for F7ext = 0 is shown in 
Fig. and allows for the definition of the critical temper¬ 
ature. We chose the temperature scale such that the orig¬ 
inal experimental temperature is 1 = T* = 298 K, i. e., 
25 °C. The phase transition then occurs at Tc = 0.86, 
and the results for the magnetization show that the phase 
at T < Tc is ferromagnetic, the phase at T > Tc para¬ 
magnetic. Returning to the binary mixture picture, the 
transition “paramagnetic-ferromagnetic” corresponds to 
“clustered-demixed”. 

To estimate the influence of the cutoff radius i^cut? 
we performed IMC simulations for image Gl#l using 
the values Rent = 12,16,..., 24 in addition to the case 
^cut = 22 discussed above. The resulting potentials 
have the same overall shape, deviations at short distances 
(r < 10 px) are small and hardly visible in a plot. The 
most significant deviations occur near the respective cut¬ 
offs, where we attribute the potentials’ fluctuations to 
statistical errors anyhow. Thus the exact value of Rent in 
the IMC algorithm is not crucial. However, although the 
potentials differ only slightly, the Tc values obtained from 
direct MC simulations scatter somewhat but with no sys¬ 
tematic trend. They lie in the range Tc G [0.81,0.86], 
corresponding to a deviation of 6%Tc. 

We further analyzed images of double- and triple-layer 
graphenes, which are shown in Appendix The results 
are qualitatively similar, regarding the potentials them¬ 
selves (Appendix [E|) as well as the simulation results: all 
transitions occured in the range Tc G [0.85,0.88] and no 
systematic dependence on the number of graphene layers 
could be observed, see Appendix \F\ Fig. 

For the different images and different i^cut, Tc G 
[0.81,0.88] was found. Converting this to kelvins, the 
transitions are predicted to occur at Tc G [241, 262] K, 


N 
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i. e., [—32, —11]°C. We note that cn of a monomolecular 
fluid film cannot be measured on the background of the 
large heat capacity of the substrate, but is a very sen¬ 
sitive numerical indicator of the transition, which is in 
turn observable in the experiment. Our estimate for 
gives a guidance for further experimental work. 

Images of the spin states and plots of correlation func¬ 
tions show that the clusters grow quickly with tempera¬ 
ture decreasing from T = 1. At temperatures increasing 
from T = 1, the cluster sizes decrease but the system 
remains in a clustered state at least up to T = 2 and pre¬ 
sumably even further. Nothing special happens at the 
experimental temperature T* = 1. 

Let us summarize our findings. We used an inverse 
Monte Carlo approach to reconstruct the effective in¬ 
teractions within a monomolecularly thin fluid film of 
a water-ethanol mixture in a slit-pore between mica and 
graphene surfaces. The shape of this potential qualita¬ 
tively confirms our qualitative argumentation in [2]. Di¬ 
rect Monte Carlo simulations with this potential allow 
for obtaining information which is not immediately ac¬ 
cessible from the experiment, for example, the critical 
temperature of the demixing transition. We considered 
two-dimensional systems, but the method can be applied 
to systems of any dimension. 

Appendix A: Binary Mixture and Lattice Spin 
Models 

Let us consider the lattice model of a binary mixture, 
in which each lattice site can be occupied by one of the 
species, A or 5, i. e., the variable characterizing the occu¬ 
pation of the site i, qi is either A or B. The interactions 
between the different species depends on the distance and 
therefore is given by the set of three distance-dependent 
interaction energies, VAA(ni), and 

where Vij denotes the distance between the correspond¬ 
ing lattice sites. 

The energy of the system with given number of A and 
B sites is then 

B = ^ [VAA{rij)6q^^A^q.^A + VBB{rij)6q.^B^qj ,B 

ijeG 

+ VAB{rij){Sq.^A^q.^B + ^qi,B^qj,A)] • (Al) 

The Kronecker deltas can be replaced by functions of 
“spin”-variables equal to -M if the site is occupied by 
A or —1 if it is occupied by B. Then hq.^A = (1 + cr^)/2 
and hq.^B = (1 “ Substituting Kronecker symbols 

for these expressions and collecting similar terms we get: 

[VAA{rij)+ VBBirij)+ 2VAB{rij)] 

ijeg 

+ ^ E E i^AAirA - VBBirA] (A2) 

ieQ jeG 

+ i E i^AAiVij) + VBBiVij) - 2VAB{rij)] (Jidj. 

ijeG 
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FIG. 5: Depiction of the concept of radially symmetric lattice 
potentials: Lattice plot of the reconstructed potential in Fig.|^ 
(Gl#l in Fig. [9). Golor plot with linear grayscale coding 
(white = high value). 

The first term is an offset which does not depend on the 
configuration and may be neglected (cancels out in the 
entropy S). Because of the site-independent interaction 
potential, the sum over j in the second term does not 
depend on i. This term, after summing up over all lat¬ 
tice sites j (for large lattices and summable potentials) 
gives a constant multiplied by the sum over cr^. The cor¬ 
responding constant is later incorporated in iLext- 

= W - VAAirA] ■ (A3) 

jeG 

The third term contains the effective interaction poten¬ 
tial, similar in spirit a Flory-Huggins parameter in the 
short-range interaction case, 

W{rA = ^ [VAAirA + VBBirA - ^VAB(rA] • (A4) 

In open systems an additional term, accounting for the 
energy difference between the molecules in the gas phase 
and in the pore, has to be taken into account. This term 
contributes to the chemical potential, is linear in the con¬ 
centrations and is incorporated in iLext- Thus, the pa¬ 
rameters of the model are a scalar iLext and the effective 
potential IT(r), for discrete values of r corresponding to 
the distances between the lattice sites. While we use a 
lattice spin model in our investigations, the discussion 
above allow for converting any quantities of interest into 
the corresponding one for a binary mixture. 

Appendix B: Details of Implementation 

The long-ranged nature of the interaction potentials 
in the system makes the direct MC simulation computa¬ 
tionally demanding. In order to keep the complexity of 
the problem acceptably low already at the stage of direct 
MC, we introduced a cutoff radius that limits the (other¬ 
wise possibly infinite) range of the pair potential W. We 
denote this cutoff radius by Rent • Formally, the potential 
is Wr = 0 for distances r > Rent- Numerically, we do not 
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sum over these distances. The values of Rent used in our 
simulations are given in Tab. [Ilj 

In order for equation to be solvable, the covariance- 
matrix of the spin-spin correlations must have full rank. 
Since it is the Hessian matrix of the entropy S' as a func¬ 
tion of IT, this should be the case near a local minimum, 
where the matrix is positive definite. However, if we let 
d be the dimension of the system of equations (§, i.e., 
the number of independent potential values to be recon¬ 
structed, then the statistical estimates for the covariance 
matrix (obtained from direct MC) may still be singular 
if they are obtained from less than d + 1 measurements. 
Hence one has to take more than d measurements in each 
IMC iteration. 

As discussed above, the complexity of the algorithm 
increases with the size of the system of equations §, 
i.e., with the number of IT values that are to be recon¬ 
structed. For this reason, one should keep the dimension 
of the algebraic system low, in order to have feasible com¬ 
puting times. Therefore we consider radially symmet¬ 
ric potentials. Of course, exact radial symmetry is not 
possible on a square lattice. And the radial symmetry 
is a postulate rather than an experimental observation. 
Nevertheless we introduce radial binning of the distances 
and consider only the radial distance instead of two- 
dimensional vectors. The input data are correlation func¬ 
tions as shown in Fig. [^with binned distances. The bins 
are chosen such that all nearest neighbors (r = 1), next- 
nearest neighbors (r = v^), axial next-nearest neighbors 
(r = 2), and so on (r = v^, ^/S, 3,...) fall into an indi¬ 
vidual bin, respectively. Hence a geometrical meaning is 
conserved as far as possible. At larger distances bins are 
merged if they are close to each other (our criterion: dif¬ 
ference < 0.4). An example of such a potential is shown 
in Fig. 

Since the individual bins have different volumes, the 
governing equations have to be adapted slightly. With 6^ 
we refer to the volume of bin p, i. e., the number of pairs 
of sites whose distance is inside bin p. The corresponding 
correlation function is 


C'p(W) ^ 

^ ideQ\rij in bin p 

The system of linear equations has to be modified 
accordingly: Ax = y with 


converged. But since the input correlation C is computed 
from a single SFM image 




(B5) 


no time-series averaging is possible as in a MC simula¬ 
tion, that would allow to obtain a statistical error from 
fluctuations. 

Hence, we need to estimate the statistical error from 
the single image we have at hand. We denote the bin 
averages in equation (Bl) by Cp{{cF}) = {{diCFj))r^jr^p> 
As an estimator for the statistical error of this average, 
we use the sample variance 


AC. 


1 2 


1 


2 


i—T (' - ■ 


(B6) 


The iterations were stopped when either the prescribed 
value of the relative error 




^ {{c,)\w-c,y 


Nb ^ 

p=l 


AC. 


1 2 


(B7) 


was achieved ('d < 1 was used), or after the maximal 
predefined number of iterations, if the convergence was 
too slow. Cases where'd < 3 could not be achieved were 
considered to not have converged. The a-posteriori values 
of achieved relative errors are given in Tab. [TT| 


Extension of the Algorithm 


As already mentioned, the external magnetic field i^ext 
can be included in the algorithm in a fully analogue way 
as the pair interaction IT. The derivative of the Hamil¬ 
tonian R with respect to i^ext is (cf. equation 0 ) 


affext 




(B8) 


To determine also i^ext, one must include the observable 
Ms as well as its variance and covariance with all the 
Cp i nt o th e simulation. The linear system of equations 
(B2)-(B4) has to be extended by 


^pq 

P Cov(6p Cp, Cg) ^ , 

(B2) 

Xp 

= 

(B3) 

Up 

= bp {Cp)\y^ — Cp 

(B4) 


V p, g G {!,..., Nb}, where Nb is the number of bins. 

To set up a suitable convergence criterion, we compute 
the weighted mean square deviation of {Cp)\w from Cp, 
where the weights are the statistical errors of the input 
data at the corresponding distance p. This deviation has 
to be of the order of unity for the algorithm to have 


ATb + I 
^Vb + 1 Vb + 1 
XNb + 1 
VNe + I 


(^Nb+Ip — P Cov{hp Cpj , 

/3 (B9) 

-(5Fext, (BIO) 

(Bll) 


V p G {!,..., Nb}, where Var(A) = Cov(A, A) and M is 
the magnetization from the SFM image. The algorithm 
will then also provide corrections in each iteration, 

when the system of size (A'^ + l) x (A'^ + l) is solved. Like 
for the SW, a damping factor A has been used with ^iLext 
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and the relative square deviation of the mag neti zation 
has been included in the convergence criterion (B7). One 
should note, that the corrections are obtained not 

only from the simulated magnetization, but also from the 
correlations at all distances. 

The empirical finding with this extended algorithm is, 
however, that the numerical stability and the speed of 
convergence may reduce notably when i^ext is allowed 
to vary (“float”). A fixed value of i^ext was thus used 
whenever possible. For images with specific magnetiza¬ 
tions m near 0 (white-pixel concentrations Cwhite near 
50%), |i7ext| is small and we obtained suitable values for 
it by trial and error. This was the case for 4 SFM im¬ 
ages. In the opposite case (2 SFM images), where m is 
not small, we could not obtain (guess) suitable values 
for i^ext: because the system’s behavior changes heavily 
as the interaction potential is iterated. In these cases, 
i^ext was kept floating in the IMC procedure. This is 
illustrated by the values given in Tab. |n| 


Appendix D: SFM Images 


TABLE I: List of the six experimental SFM images analyzed 
in this work, crfiiter is the width of the Gaussian convolution 
filter that was applied to the (flattened) raw image. Cwhite 
and m are the white pixel percentage and the magnetization 
per volume for a spin system, respectively. Gl^l is the image 
discussed in the main text. 


* Image G3#l was not filtered, only flattened and reduced 
to Ibit. It was scaled down from 512 pixel beforehand to 
simplify the inverse MG procedure. 


name 

layers 

width/px 

•^filter/P^ 

width/nm 

Cwhite /% 

m/% 

Gl#l 

1 

512 

1.0 

500 

48.6 

-2.9 

Gl#2 

1 

512 

1.0 

533 

67.0 

+34.0 

G2#l 

2 

512 

1.0 

1000 

38.4 

-23.2 

G2#2 

2 

512 

1.0 

1000 

51.4 

+2.8 

G2#3 

2 

512 

1.0 

1000 

51.4 

+2.7 

G3#l 

3 

256* 

0.0 

1000 

53.0 

+5.9 


Appendix C: Test case: Ising model with 
dipole-dipole interaction 


We tested the implementation of the IMC algorithm 
on a two-dimensional ferromagnetic Ising model with ad¬ 
ditional long-ranged repulsive 1/r^ interaction, as dis¬ 
cussed before in the context of ultrathin magnetic films 
[4] which exhibit a clustered state at intermediate tem¬ 
peratures. The Hamiltonian of this system is given by 

'Htest({cr}) = -^ y] CTiCTj 

+5 E (Cl) 


The first term is a nearest-neighbor interaction (Ising 
model, ferromagnetic for J > 0), the second term an 
out-of-plane dipole-dipole interaction (repulsive for par¬ 
allel spins), the third term is due to an external magnetic 
field. This Hamiltonian corresponds to a pair potential 



J 

2 


1 for r = 1 
0 else 


(C2) 


with r = ||r||. 

The results in Fig. show that the IMC procedure 
is able to reconstruct the underlying pair-potential from 
CFs obtained from a direct simulation. For larger dis¬ 
tances and hence small values, the reconstruction is lim¬ 
ited by noise. If the number of MC steps in the direct 
simulation of the input CF is increased, this noise level 
decreases due to the improved statistical precision of the 
input data. Since (unlike for experimental data) the sta¬ 
tistical error of the input CF can be precisely measured 
in the direct MC simulation, 'd < 1 is reached by the 
IMC, cf. Eq. (1^). 


The experimental SFM images that were analyzed in 
this work are briefly described in Tab. |Tj The images 
themselves are shown in Fig. 

Image G2^3 shows a dark, slightly curled line from 
bottom middle to top left. This is a depression in the 
graphene membrane of unclear origin. Possible causes: 

• a defect in the mica crystal, 

• a defect in the graphene sheet, 

• a wrinkle (downward!) in the graphene possibly 
caused by compressive stress. 

In any case it noticeably affects the clustering patterns. 

Gray pixels in the images indicate void areas, that were 
not included in the correlation functions. 


Appendix E: IMC Results 

In Fig. we compare the original SFM image with a 
spin configuration from a MC simulation with the recon¬ 
structed potential. Even though any quantitative judg¬ 
ment should involve the correlation function, this depicts 
the similarity between experimental and reconstructed 
system. 

The reconstructed potentials for all six SFM images are 
shown in Fig.[^ The parameters used for the IMC proce¬ 
dure are listed in Tab. [Til along with the computing times 
and the deviations from the experimental data that was 
reached. Remarkable features that all potentials have in 
common: short-ranged attractive, then somewhat longer- 
ranged repulsive and apparently fluctuating around 0 for 
larger distances. The potentials from Gl#l-G2#2 show 
a very similar pattern: the first two values are negative, 
where Wi ~ 1.8 then four positive values in a char¬ 
acteristic zigzag pattern with the most repulsive value at 
the third shortest distance r = 2 then apparently unsys¬ 
tematic fluctuations around LF = 0. Deviations occur for 
G 27 ^ 3 , which has a dark line across the image (cf. Fig.[^ 
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FIG. 6: (Color online) Reconstructed potential for the “Ising + dipole-dipole” test case Eq. (Cl) with J = 10, ^ = 1, 
= 0.05, T = 10. The input CF was generated in 10000 MC steps, with measurements every 13 steps. L = 64 and 
Rent — 16 were used in both, the direct and the inverse MC algorithm. A = 0.10. A random initial potential of small amplitude 
was “preconditioned” by local potential corrections (6996 MC steps) [TS] . The correlation Cr converged after 19 inverse MC 
iterations = 50171), yielding a deviation of= 0.99, cf. Tab. \n\ The reconstructed external field is i7ext = 0.0477. 

The double-logarithmic plot in the lower panel reveals that the ijr^ behavior can be reconstructed. Negative values in Wr 
occur at r = 1 and r > 9.87 (not shown in the log-log plot). 


and the triple-layer image G3#l, which has the largest 
structures of all images. 

At this point, we would once again like to point out the 
geometrical meaning of distances on the lattice. Fig. 
shows the reconstructed potential from Gl#l on the lat¬ 
tice for short distances. One can clearly see a square¬ 
shaped potential barrier. This certainly has a geometri¬ 
cal meaning connected with the lattice geometry and the 
implementation of certain concepts on a lattice, e. g., an 
elastic line-tension that we expect to be present due to 
the elastic graphene deformation at the cluster bound¬ 
aries. 

Furthermore, the results are effective potentials. One 
should thus not interpret such a result as: “The mi¬ 
croscopic, inter molecular potential is attractive at a dis¬ 
tance of \/2 pixel = 1.38 nm and repulsive at a distance 
of 2 pixel = 1.95 nm.” 


Possible Error Sources 

There remains some uncertainty about the scope of 
application of the reconstructed potentials. First of all, 
the SFM images are prone to noise and thermal drift. 
Hence we corrected for the drift, applied a Gaussian fil¬ 
ter and reduced the color depth. This manipulation to 
the raw image could introduce an error. However, we an¬ 
alyzed the influence on the correlation function [2] with 
the finding that it is quite stable to the image manipula¬ 
tion. 

Secondly, the statistical basis obtained from one im¬ 
age (512 X 512 pixel) could be too poor to extract even 
the main characteristics of the interaction potential. In 
addition, our estimate of the input correlation function’s 


error is very coarse. 


Appendix F: MC Study of the Resulting System 

In direct MC simulations using the reconstructed po¬ 
tential, we examined the three observables specific heat 
capacity, (specific) magnetization, and magnetic suscep¬ 
tibility 


Cn 


m 


Xn 


1 rdE\ _ {n^) - {uf 

N \dT) ~ NT^ 

M I ^ 

~N ~ 

ieQ 

1 r dM \ _ (M^) - (m)2 

N ~ ^ 


(FI) 

(F2) 

(F3) 


We used the potential Gl#l to obtain all the data 
shown here. However, simulations with G2#l and G3#l 
yielded similar results. 


1. Variable Field Hext 


We fixed T = 1 and ramped F7ext up and down. 
No hysteresis could be detected in the range F7ext ^ 
[0.00,0.15]. The magnetization m sensitively increases 
with ii^ext and shows saturation already at values below 
^ext = 0.1. Correspondingly, the susceptibility xn is 
high. The external field depends on the chemical po¬ 
tential in the binary mixture picture. Hence a differ¬ 
ence between Vaa and Vbb (iu the sense of Eq. (A3)) by 
less than 0.2 /cpT leads to a system dominated by one of 
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FIG. 7: The six experimental SFM images analyzed in this work. Row-wise from top left to bottom right: Gl#l, Gl#2, G2#l, 
G2#2, G2#3, G3#l. Flattened, filtered and reduced to 1 bit. Void pixels are shown in gray. Gl#l is the image discussed in 
the main text. 
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Detail of SFM image Gl#l, 128pixel or 125 nm wide. Snapshot of MC simulation with the corresponding re¬ 
constructed potential and i/ext = —0.001 on a 128 x 128 
lattice. 

FIG. 8: Gomparison of experimental image Gl#l and simulation with the reconstructed potential for this image. 
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FIG. 9: (Golor online) Reconstructed potentials from SFM images: single- (Gl#n), double- (G2^n) and triple-layer (G3#l) 
graphene. Data is gradually shifted upward to improve readability. The corresponding zeros are indicated by dashed horizontal 
lines. Note that the r-axis is in pixels and that images have different resolutions, cf. Tab. [I] and Fig. 


TABLE II: Parameters (middle part) and convergence data (right part) for the inverse Monte Garlo algorithm, niter is the 
number of IMG iterations after which convergence oc curr ed, is the corresponding number of global MG updates, is the 
relative correlation/magnetization deviation (cf. Eq. (B7)) at the next measurement, tcpu is the GPU time on one thread of 
an Intel® Gore^^ 17-3770 GPU @ 3.40GHz. 

The iLext values marked with an asterisk (*) were fixed, the others are results from the algorithm. The specific magnetization 
m (from Tab.|I]) in the second column is shown once more to illustrate where iLext can be fixed (|m| small) and where not {\m\ 
large). 


image 

m/% 

L 

Rcut 

77ext 

A 

niter 

U/10" 


tcpu 

Gl#l 

-2.9 

64 

22 

-0.001* 

0.50 

14 

7.9 

0.98 

12h 

Gl#2 

+34.0 

80 

20 

+0.0056 

0.10 

14 

3.7 

1.96 

8h 

G2#l 

-23.2 

80 

24 

-0.0035 

0.15 

23 

19.5 

1.99 

56 h 

G2#2 

+2.8 

80 

24 

+0.001* 

0.15 

14 

6.3 

2.00 

18h 

G2#3 

+2.7 

80 

20 

+0.000* 

0.10 

27 

11.8 

1.05 

24 h 

G3#l 

+5.9 

89 

22 

+0.001* 

0.50 

30 

5.7 

2.85 

18h 
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temperature T temperature T 

FIG. 10: (Color online) MC simulations with reconstructed potential (Fig.[^: Observables m (left panel) and xn (right panel) 
in dependence on T with fixed i/ext = 0.0. Measured with decreasing T in four independent series (different markers) on a 
64 X 64 lattice, 10^ MC updates per measurement (of which 10% for thermalization; T G [0.81, 0.87] was resolved with 10^ MC 
updates per measurement). 


the components. Back in the spin picture, we are obvi¬ 
ously dealing with a paramagnetic phase here (i. e. near 
T = 1). 

2. Variable Temperature T 

The temperature dependence of c^r for i^ext = 0 is 
shown in Fig. an d we add here the dependence of m 
and Xat in FigT]lQ[ Snapshots of spin configurations at 
different temperatures are shown in Fig. El One can 
clearly see the phase transition happening between T = 
0.84 and T = 0.87, where the black-white symmetry is 
broken which means spontaneous magnetization. 


3. Results for Thicker Graphenes 

Apart from the potential for image Gl#l (single-layer 
graphene; m ^ 0, i. e., black/white symmetric) discussed 
in the main text, we also performed MC simulations with 
the potentials for images Gl#2 (single-layer, m > 0), 
G 27^1 (double-layer, m < 0) and G3#l (triple-layer, 
m « 0). In all cases the external field was set to i^ext = 0 
and the temperature was decreased from T = 1 to ob¬ 
serve the critical point. Results for the thermodynamic 
observables are shown in Fig. The critical temper¬ 
atures (Te = 0.88 for Gl#2, G [0.86,0.88] for G2#l 
and Tc = 0.85 for G3#l) show no systematic dependence 
on the number of graphene layers. 
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FIG. 11: Snapshots of spin states from MC simulations with the reconstructed potential for image Gl#l. Selected temperatures 
T; 64 X 64 lattice, cf. Fig.[^ 





temperature T 


FIG. 12: (Golor online) MG simulations for the potentials (cf. Fig. reconstructed from images Gl#2, G2#l and G3#l 
(single-, double- and triple-layer, respectively): Observables cn (left panel), m (middle) and xn (right panel) in dependence 
on T with fixed i7ext = 0.0. Measured with decreasing T. The absolute value of the averaged magnetization |(m)| is plotted 
instead of the magnetization (m), because its sign is arbitrary in the absence of an external held. 

The critical temperatures are Tc = 0.88 for Gl#2, Tc G [0.86,0.88] for G2#l and Tc = 0.85 for G3#l, respectively. 



















































































